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O ■ 1 INTRODUCTION 

In the last decades, it has become clear that energetic feed- 
back processes are the key ingredients in shaping the star 
formation history of the Universe and regulating the for- 
mation and evolution of galaxies. In particular, supernovae 
(SNe) and black holes (BHs) are considered the principal 
sources of such feedback (see Silk 2007, for a review). 

Cosmological simulations have been one of the most 
powerful tool for investigating the formation and evolu- 
tion of galaxies (see review from Bertschinger 1998; Springel 
2010). However, current cosmological simulations lack both 
the resolution and the physics needed to model small scale 
feedback phenomena like the explosion of a single SN or the 
accretion into a BH. Several numerical, sub-grid recipes have 
been employed to mimic the large scale effect of feedback by 
injecting the energy either in thermal or kinetic form. 

Thermal feedback has the advantage of having an 
isotropic effect on the closest surrounding medium. How- 
ever, in the presence of radiative cooling and poor res- 



ABSTRACT 

We perform simulations of feedback from supernovae and black holes with smoothed 
particle hydrodynamics (SPH). Such strong perturbations are inaccurately handled 
with standard time integration schemes, leading to poor energy conservation, a prob- 
lem that is commonly overlooked. We show for the first time that, in the absence 
of radiative cooling, concordance of thermal and kinetic feedback are achieved when 
using an accurate time integration. In order to preserve the concordance of feedback 
methods when using a more efficient time integration scheme - as for instance the 
hierarchical time-step scheme - we implement a modified version of the time-step lim- 
iter proposed by Saitoh & Makino (2009). We apply the limiter to general test cases, 
and first show that this scheme violates energy conservation up to almost four orders 
of magnitude when energy is injected at random times. To tackle this issue, we find 
necessary, not only to ensure a fast information propagation, but also to enforce a 
prompt response of the system to the energy perturbation. The method proposed here 
to handle strong feedback events enables us to achieve energy conservation at percent 
level in all tests, even if all the available energy is injected into only one particle. We 
argue that concordance of feedback methods can be achieved in numerical simulations 
only if the time integration scheme preserve a high energy conservation level. 
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olution, a large part of the energy is radiated away be- 
fore any hydrodynamical response of the medium. Sev- 
eral studies concentrated their efforts in numerically solv- 
ing the over-cooling problem and the consequent dissi- 
pation of thermal feedback energy (e.g. Gerritsen 1997; 
Mori et al. 1997; Thacker & Couchman 2000; Kay et al. 
2002; Sommer-Larsen, Gotz, & Portinari 2003; Brook et al. 
2004; Stinson et al. 2006; Booth & Schaye 2009). In order 
to avoid the over-cooling problem, other studies favoured 
the kinetic feedback approach, in which all or part of the 
energy is given to the surrounding medium as momen- 
tum (e.g. Navarro & White 1993; Mihos & Hernquist 1994 
Kawata 2001; Kay et al. 2002; Springel & Hernquist 2003 
Oppenheimer & Dave 2006; Dalla Vecchia & Schaye 2008 
Dubois & Teyssier 2008). However, over-cooling can still 
play an important role in quickly dissipating the heat pro- 
duced by the outflowing gas shocks, and underestimating 
the increase of the temperature of the medium. 

Until now, despite the large variety of implementation 
of feedback models, very little work has been done to com- 
pare the two feedback approaches in the context of galaxy 
formation. Kay et al. (2002) compared the results of cosmo- 
logical simulations with either thermal or kinetic SN feed- 
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back. They showed that even when artificially preventing 
cooling of heated particles the results of the two schemes 
are considerably different. Nevertheless, before tackling the 
over-cooling problem, one should first make sure that the 
input energy is accurately conserved. 

In a recent study of the modelling of SN feedback in 
SPH simulations, Saitoh & Makino (2009) (hereafter SM09, 
see also Merlin et al. 2010) noted that strong perturbations 
of the internal energy of gas particles lead to the violation 
of energy conservation when an individual particle time-step 
scheme is used (see Fig. 1 for an illustration). That can 
be summarised as follows. If a source of energy alters the 
internal energy of a (active) gas particle, the particle will 
eventually decrease its time-step according to the Courant 
criterion. However, if its neighbouring particles are inactive,^ 
they will not react immediately to the change of the thermal 
state of the region, and the integration accuracy may be 
strongly compromised. The lack of a prompt response to the 
nearby energy injection can lead to effects like inter-particle 
crossing (due to particles missing the deceleration given by 
viscous forces) and non-conservation of energy. 

SM09 proved that, when using an appropriate time- step 
limiter, accurate energy conservation is achieved for strong 
internal energy perturbations. They proposed a scheme in 
which the ratio of the time-steps (longer over shorter) of 
neighbouring gas particles cannot be larger than a given 
factor /step- A value of /step = 4 is enough to ensure energy 
conservation to similar level of accuracy given by a global 
time-step integration scheme (Merlin et al. (2010) suggested 
the values /step = [4, 8]). 

In this work and for the first time, we investigate un- 
der what conditions thermal and kinetic feedback schemes 
lead to qualitatively identical results. We first show that the 
feedback prescriptions are equivalent if the integration of the 
hydrodynamical equations is done over global, system time- 
steps. In order to do that, we compare test simulations of 
Sedov's blast wave problem (Sedov 1959; Landau & Lifshitz 
1959) to its analytic solution, finding very good agreement 
and energy conservation to within 1%. 

With the aim of generalising the results, we investigate 
whether one can achieve similar agreement and accuracy 
with the more computationally efficient individual time-step 
integration scheme. We implement a modified version of the 
SM09 time-step limiter in the public version of GADGET-2 
(Springel 2005) consistently taking into account the time- 
step synchronisation and the underlying leapfrog integrator. 
We test the robustness and speed of the limiter with Sedov's 
explosion problem injecting either thermal or kinetic energy 
at a random time. We demonstrate that thermal and kinetic 
feedback methods are equivalent and energy is accurately 
conserved, provided that the system promptly responds to 
the energy perturbation. 

With the application of the time-step scheme to galaxy 
formation and cosmological simulations in mind, we study a 
more realistic problem: the release of energy in the presence 
of density and pressure gradients. We apply the scheme to an 
off-centre explosion in a self- gravitating gas halo, and show 
that the concordance of feedback methods and accurate en- 
ergy conservation are still achieved even for the extreme case 
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Figure 1. Sedov's explosion test for thermal (left column) and 
kinetic (right column) energy injection, when the individual time- 
step integration is used. Results are given in the natural system 
of units soon after the explosion at time t = 0.004. We show in 
the upper panels the projected gas density in a slice of thickness 
Az = 0.1 centred on the box origin. The dashed circle corresponds 
to the position of the spherical shock-front given by the similarity 
solution in Eq. 2 at that time. The white dots mark the position 
of all particles that initially received the energy. In the lower 
panels, the average radial density profiles are compared to the 
profiles given by Sedov's analytic solution (red lines). Given the 
large energy jump introduced by the explosion, it is clear from 
these plots that none of the feedback scheme are able to describe 
properly Sedov's test. 



where all available energy (either thermal or kinetic) is in- 
jected into only one particle. 

In standard cosmological simulation practice the new 
integration approach presented in this paper could, in princi- 
ple, require more simulation steps, but will mainly populate 
the intermediate time-bin levels with more particles. There- 
fore, we could expect a slight increase of the computational 
time for galaxy formation simulations, a price to pay for a 
better handling of feedback events. However, we still want 
to emphasise that, in idealised simulations like the ones pre- 
sented in this work, the proposed scheme enables much faster 
integration than the standard, individual time- step scheme 
since it preserve a high energy conservation accuracy. 

The paper is organised as follows. We first show in Sec. 2 
that concordance of feedback methods is achieved with a 
global time-step scheme. We then apply our individual time- 
step algorithm to Sedov's blast wave test and to the off- 
centre explosion in a self- gravitating gas halo in Sec. 3. We 
conclude in Sec. 4. 

In all the simulations presented in this paper, we use 
the so-called natural system of units in which the units of 
velocity, length and mass are unity. 
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Name 


Scheme 




a 


Energy (%) 


Time (h) 




Thermal 


Kinetic 


Thermal 


Kinetic 


SEDOV-O 


Olobal 


1 


2 


0.81 


0.78 


3.71 


2.50 


SEDOV-G-£;01 


Global 


0.1 


2 


0.83 


0.82 


1.97 


1.46 


SEDOV-G-£;001 


Global 


0.01 


2 


0.84 


0.74 


1.09 


0.96 


SEDOV-G-a 1 


Global 


1 


1 


1.99 


2.30 


3.87 


2.73 


SEDOV-G-a4 


Global 


1 


4 


0.73 


0.61 


3.58 


2.38 



Table 1. Numerical parameters and energy conservation estimate and computational time for Sedov's blast wave tests with a global 
time-step scheme. From left to right: names of the runs specifying the integration scheme and the change of numerical parameters, time 
integration scheme; injected energy, E^-^ artificial viscosity, a; energy conservation estimates in percent (from Eq. 3); and computational 
time in hours. Statistics are given at the end of the runs, at time t = 0.1. The wall-clock time is given for runs on 8 cores. The reference 
simulation is highlighted in bold. 



2 CONCORDANCE UNDER ACCURATE 
TIME INTEGRATION 

In this section we demonstrate the concordance of thermal 
and kinetic feedback schemes. We run several simulation 
tests of Sedov's explosion problem varying the numerical 
parameters to show the robustness of our findings. Here, 
we assume the most favourable case in term of integration 
accuracy. The integration is done over global time-steps in 
order to concentrate only on how and why concordance is 
achieved. 

We recall that Sedov's problem with thermal energy 
input has already been numerically solved with SPH (e.g. 
Springel & Hernquist 2002; Tasker et al. 2008), reaching the 
desired integration accuracy by limiting the particle maxi- 
mum time- step or evolving the system on small, constant 
time-steps. 

2.1 Simulation setup 

We generate glass-like initial conditions (White 1996) with 
N — 128"^ gas particles by evolving a random distribution of 
particles with inverted gravitational force sign. We further 
relax the glass particle distribution by running it with hydro 
forces only, in order to smooth the pressure gradients that 
may remain due to the small density fluctuations. In the 
natural system of units the size of the box is L = 1. In order 
to obtain a uniform background density, po = 1, the gas 
particle mass is set to rrig — 1/N. We assign to each particle 
the internal energy uq = 10~^ so that the late evolution of 
the explosion is not affected by the energy of the medium 
that the expanding shell accumulates with time. 

The total energy is injected in a small region in the 
centre of the volume either in thermal or kinetic form. In 
the simulations presented in this section, we heat/kick the 
^[h,k] = 32 particles within a sphere of radius ro (we refer 
to heated and kicked particles with the subscripts 'h' and 
'k', respectively). The energy is distributed using the SPH 
kernel weight normalised to unity, and with kernel size ro 
Particle i at distance from the centre of injection receives 
the energy per unit mass, 

_ w{ri,ro) E^ 



where rrig is the mass of the particle, and w{ri,ro) is the 
kernel weight at distance r^. To conserve energy, the ker- 
nel weight is normalised to the sum of the weights over all 
heated/kicked particles. 

If the particle is heated, its internal energy increases by 
Ui. If the particle is kicked, its velocity is increased by Vi — 
\j2ui^ and the change of kinetic energy in this case equals the 
change of thermal energy in the other. The momentum kick 
is radial, along the direction from the centre of the explosion 
to the particle. 

For a given number of heated/kicked particles, because 
the injected energy per unit mass is u^, (x E^/rrig oc N ^ 
the mass resolution defines the energy jump. Our choice of 
initial conditions leads to an energy jump of the order of 
10^ with respect to the medium initial energy and for the 
reference value E^ — 1. 

For given values of E^ and po , and assuming monatomic 
gas with a ratio of specific heats 7 = 5/3, the analytic solu- 
tion for the time evolution of the shock-front radius is given 
by (Sedov 1959): 

1/5 



R{t) ^ 1.1527 



,2/5 



(2) 



(1) 



Regarding the time integration, we employ in this sec- 
tion global, adaptive time-steps. At each simulation step, 
the whole system is advanced in time by integrating over 
the minimum time-step of all particles. To define the time- 
steps, we use the standard Courant parameter, C — 0.15, 
and the value of the time accuracy coefficient 77 = 0.0025 in 
equations Bl and B3 presented in Appendix B. Regarding 
the accuracy parameter 77, we performed accuracy and effi- 
ciency tests with different values, and chose 77 = 0.0025 as 
fiducial parameter for this study. For a comparison of those 
tests, the reader should refer to Sec. 3.2 and 3.3. Finally, 
in order to define the SPH kernel, we set the number of 
neighbours to 32 di 2 particles. 

The list of runs is given in Table 1. In the reference 
model (highlighted in bold) we adopt the values E^ — 1 
and a — 2 for the injected energy and the artificial viscosity 
coefficient, respectively. We also run a sample of simulations 
varying E^ and a. Beside the numerical parameters of the 
tests, we also list in Table 1 a measure of the conservation 
of the injected energy at the end of the runs: 

A^. _ {E+ - Eo) - 



(3) 
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Figure 2. Sedov's explosion test for the reference simulations, 
where the system is evolved on global time-steps: energy is in- 
jected either in thermal (left column) or kinetic (right column) 
form. Results are given at t = 0.02 after the explosion. We see 
that both methods agree remarkably well in reproducing the sim- 
ilarity solution for both the location and the width of the blast 
shell. 



where E-^ = + is the time evolution of the sum of the 
thermal and kinetic energies in the simulated volume and 
^0 is the initial, total energy of the system (excluding the 
input energy) . We also list the total simulation running time 
on 8 cores. 



2.2 Simulation results 



We first describe the reference model by showing in Fig. 2 
the blast wave behaviour at time t = 0.02. In the upper row, 
the projected gas density of a slice of thickness Az = 0.1 cen- 
tred on the box origin is plotted for the thermal and kinetic 
feedback runs. The white dots mark the position of all parti- 
cles that initially received the energy. The dashed circle cor- 
responds to the position of the spherical shock-front given by 
the analytic solution in Eq. 2 at that time. In the lower row, 
the simulated, average radial density profiles (black lines) 
are compared against Sedov's similarity solution (red lines). 

The agreement with the analytic solution is very good 
for the two feedback schemes, for both the location and the 
width of the blast shell. As shown by other authors, the lack 
of resolution and the SPH smoothing do give a density con- 
trast at the shell radius that is smaller than the theoretical 
prediction (see for instance a comparative study for AMR 
and SPH codes by Tasker et al. 2008). Nevertheless, we see 
on the radial profiles that our resolution is sufficient enough 
to describe correctly the asymmetry of the shock front in 
both cases. 

As expected from a global time-step integration, the 
computation time is large, ranging from ~ 3 to ^ 4 hours. 



With both feedback schemes, the accuracy of energy con- 
servation arises because at each simulation step the entire 
system is integrated. Therefore, all particles are aware of the 
hydrodynamical state of their neighbours. We will develop 
this argument in more detail in Sec. 3. Here, we focus on 
how and why there is agreement between the two schemes. 

It is somehow expected that both methods give the 
same results. Indeed, Sedov's initial conditions are the to- 
tal energy of the explosion and the medium density. The 
only requirement is that a large amount of energy is instan- 
taneously injected in a small volume,^ but its form is not 
specified (Landau & Lifshitz 1959). One can input either all 
thermal, all kinetic or a combination of both forms and ob- 
tain the same similarity solution at any given time. Since 
the hydrodynamics conservation laws are used to derive the 
solution, a property of the similarity solution is that the 
fractions of thermal and kinetic energies of the blast are 
constant in time. 

However, in the numerical integration, a finite time is 
required to convert one form of energy to the other and reach 
the energy budget given by the analytic solution (see illus- 
tration in Fig. 4). As long as momentum can be converted 
into thermal energy by physical processes like e.g. viscosity,"^ 
the numerical integration of the conservation laws should 
reach the similarity solution. 

We show in Fig. 3 the time evolution of the energy con- 
servation relative error (given by Eq. 3) for all test simu- 
lations, where solid and dashed lines refer to the thermal 
and kinetic energy injection methods, respectively. Lines of 
the same colour are for simulations with the same numerical 
parameters. Energy injection happens at time t = 0. 

We first concentrate on the reference simulation (black 
lines). We show in the plot that the violation of energy 
conservation happens in the early stage of the explosion 
{t < 10"^), when the energy contrast is the largest. In the 
thermal case (solid line), there is initially a jump of ^ 0.15% 
around t = 10~^ . The conversion of energy into momentum 
happens very quickly, and after t = 4 x 10""^ the evolution 
follows the kinetic case with an offset slowly decreasing. At 
later times, both curves flatten to roughly 0.8%. 

The same behaviour can be seen in the energy variation 
tests (blue and green lines), where the input energy is de- 
creased by a factor of 10 and 100, respectively. Decreasing 
the input energy slows the blast evolution and gives the time 
offset seen the plot. With constant a, varying the input en- 
ergy gives similar relative errors, providing an estimate that 
is independent of the energy value. Moreover, energy con- 
servation is achieved at a comparable level for both thermal 
and kinetic feedback. 

We show in Fig. 3 that the choice of the artificial viscos- 
ity parameter plays an important role in the input energy 
conservation. We compare three models with a = 1 (red), 
2 (black) and 4 (orange). The simulations with a = 1 keep 
an excess of momentum for longer time (see discussion of 



^ It has been proved that injecting energy in a region of the size of 
the spatial resolution gives an accurate treatment of the problem 
(e.g. Tasker et al. 2008). 

^ In SPH simulations kinetic energy is converted to ther- 
mal energy by the numerical, artificial viscosity scheme 
(Monaghan & Gingold 1983; Balsara 1995). 
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Figure 3. Time evolution of the input energy conservation for 
the Sedov's blast wave tests with global time-step integration. 
Solid lines show the results for the thermal injection of energy. 
Dashed lines represent the evolution of the kinetic feedback runs. 
The colour coding is defined by the properties of the runs. 



Fig. 4 below), which gives a larger energy conservation er- 
ror. On the other hand, simulations with higher a flatten to 
a smaller energy violation value (below 1%) at earlier times. 
We argue that the error accumulates during the conversion 
from kinetic to thermal energy in the shock front. This in- 
dicates that a higher artificial viscosity preserves a better 
accuracy in shocks. 

While the evolution of the input energy relative error 
informs us on the accuracy of the simulations, additional 
information is needed to understand the behaviour of the 
blast. The time evolution of the conversion of one form of 
energy to the other following the explosion is described in 
Fig. 4. In this figure, we show an estimate of the energy 
partition as [E- — Eo-)/{E-^ — ^o+), where E- — Et — E\, 
and Ej^ — Et -\- E\^ are the time evolution of the difference 
and the sum of the thermal and kinetic budget of the whole 
system, respectively. In order to follow only the conversion 
of the input energy E^^ we thus remove the initial (before 
injection) values Eq- and ^o+- If E^ is in thermal form 
(solid lines), the ratio is unity at injection time t = 0; if 
is in kinetic form (dashed lines) the ratio is initially —1. 

We compare the results with the energy partition ex- 
pected value given by Sedov's similarity solution. At a given 
time, we calculate the radial profiles of density, velocity and 
pressure, and integrate them numerically to obtain the to- 
tal thermal and kinetic energy within the blast radius.^ The 
integration gives 71.7% of the blast energy in thermal and 
28.3% in kinetic form (see Chevalier 1974, for the same re- 
sult). We plot the analytic ratio (dotted line) and find agree- 
ment with our results. 

Indeed, the evolution shows that energy of one form is 



^ Note that in Sedov's solution the total injected energy is within 
the blast radius. 



Figure 4. Time evolution of the energy partition for the Sedov's 
blast wave tests with global time-step integration. E- and Ej^ are 
the time evolution of the difference and the sum of the thermal 
and kinetic budget of the whole system, respectively. Quantities 
with subscript '0' correspond to the initial values, before energy 
injection. Solid lines show the results for the thermal injection of 
energy. Dashed lines represent the evolution of the kinetic feed- 
back runs. The colour coding is as in Fig. 3. 



quickly converted into the other, and that the ratio tends to 
the same energy partition limit value. However, as already 
mentioned. Fig. 4 shows clearly that a finite time is needed 
in order for the simulations to converge to the similarity so- 
lution. We see from the results of the different setups that 
this convergence time depends on both the physical prop- 
erties of the explosion and the numerical parameters of the 
simulation. This property will be studied in more detail in 
an upcoming work. We will especially investigate the effect 
of radiative cooling processes and compare the numerical 
convergence to the timescale at which the adiabatic Sedov's 
blast evolves to the radiative snowplow phase. 

Nevertheless, since decreasing the input energy makes 
the pressure gradient shallower (thermal feedback) and the 
viscosity term smaller (kinetic feedback) for a given artificial 
viscosity parameter, the offset seen in the time evolution can 
already be explained by the lower conversion rates. However, 
we see from Fig. 4 that changing the conversion rate clearly 
modifies the convergence process. With smaller a, viscous 
forces are reduced and more momentum is created in the 
thermal feedback case. Conversely, less momentum is dissi- 
pated in the kinetic case, even though this excess of momen- 
tum from both feedback methods tends, slowly with time, 
to the same expected limit. On the opposite, a higher con- 
version rate {a = 4) produces an excess of thermal energy, 
creating too few (thermal case) or converting too much (ki- 
netic case) momentum. With this higher a the convergence 
to Sedov's energy partition happens again slowly with time. 
We confirm here that the stability of the numerical scheme 
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depends on the artificial viscosity normalisation.^ We also 
demonstrate that the value of a = 2 reproduces best Se- 
dov's test through a faster numerical convergence. This is in 
agreement with the result of Monaghan (1992) for the choice 
of a when considering shock fronts. However, we cannot yet 
advise for this somewhat higher value than typically used in 
cosmological simulations since, in the artificial viscosity for- 
malism of GADGET-2, the a parameter does not differentiate 
shocks from shear flows. In the latter, a too high conversion 
rate would affect dramatically the velocity field at the con- 
tact surface. In this specific matter of the artificial viscosity 
in SPH simulations, we assess that more work is required 
before obtaining a more general criterion. 



3 CONCORDANCE UNDER EFFICIENT TIME 
INTEGRATION 

In simulations with a large dynamical range (e.g. astrophys- 
ical problems of galaxy and structure formation), the re- 
quired time scales can span several orders of magnitude. In 
state-of-the-art cosmological simulations, up to several bil- 
lions of particles are integrated over time. The use of con- 
stant and global time integration steps is then prohibitive in 
terms of computational costs. The individual time-step in- 
tegration scheme (first introduced by Aarseth 1963; Makino 
1991) allows particles to be integrated on time-steps which 
are functions of the local state. 

SM09 demonstrated that an accurate description of 
feedback processes which involve strong energy perturba- 
tions can only be achieved by ensuring fast information 
transfer when using individual time- steps. They proposed 
an innovative scheme in which the time-step length of neigh- 
bouring gas particles is constrained by a limiter. SM09 val- 
idated their time-step limiter algorithm with Sedov's blast 
wave test. They also tested the effect of the time- step lim- 
iter with a SN-like explosion in a self- gravitating halo of cold 
gas. They showed that the conservation of energy and linear 
momentum agree remarkably well with those obtained with 
a more conservative (but computationally more expensive) 
global time- step scheme. 

However, SM09 used initial conditions that included 
the explosion energy. This leads to setting the time-step of 
heated/kicked particles and limiting the time-step of their 
neighbours at the start of the simulation. Therefore, the in- 
tegration was correctly performed over the appropriate time- 
step from the beginning. In the most general case, e.g. in cos- 
mological simulations where the feedback events from SN or 
AGN activity do not affect only active particles, this would 
not happen. We specifically design our tests to consider the 
energy input after the simulation has been running for sev- 
eral steps. This prevents any time-step adjustment before 
the integration of the dynamical and hydrodynamical equa- 
tions is performed. We also state here that we do not limit 
the maximum size of the time- step. This setup is justified 
by the aim of describing the more general case of energetic 



^ This supports our claim that the error builds up mostly in the 
viscous term. 



feedback, where individual time-steps reflect only the hydro- 
dynamical state of particles. In order to study an asymmet- 
ric problem, we modify the case of the self- gravitating gas 
halo by off-setting the explosion. 

In the following sections, we briefly present our method 
to maintain a high accuracy when considering strong feed- 
back events in simulations using an individual time-step 
scheme. The reader interested in the implementation of this 
method should refer to the appendices for a detailed descrip- 
tion. The results of the two sets of test simulations are then 
presented, focusing on the conditions needed to achieve the 
concordance between feedback methods. 

3.1 Individual time-stepping for feedback 

In the context of the hierarchical time-stepping scheme pre- 
sented in Appendix A, it is important to note that if the 
energy content of a gas particle is suddenly increased, the 
particle itself will be informed as soon as it becomes ac- 
tive. Within a few active time- steps, a fraction of the en- 
ergy is efficiently converted from one form to the other 
(thermal to kinetic or vice versa). If thermal energy is in- 
jected, it will be converted into momentum through strong 
hydro accelerations due to the large pressure gradient. If 
heated particles and their neighbours do not adjust their 
time-steps, the integration would lead to very large veloci- 
ties and an artificial excess of kinetic energy. In the other 
case, the injected kinetic energy is converted into thermal 
energy through shocks. Eventually, the integration over a 
long time-step would lead to an excess of thermal energy, 
again violating energy conservation. 

In both scenarios, the violation of energy conservation 
is due to the fact that the particles do not react soon enough 
to the sudden change of their hydrodynamical state. Even 
if the particle becomes active at the proper time, if the new 
energy state is not taken into account in defining the length 
of the next step, the integration following the energy increase 
will also be done over a too long time- step, leading again to 
non-conservation of energy. In any case, it is therefore crucial 
to capture the initial stage of energy injection. 

We emphasise here the problems that one may en- 
counter when implementing feedback modules in SPH codes. 
We show in Appendix B that the signal velocity and the 
acceleration are the key information needed to define the 
time-step. In the public release of GADGET-2 both quan- 
tities are calculated during the computation of the local 
hydrodynamical forces. Consequently, any later change of 
the energetics of an active particle would not be taken into 
account during the calculation of the next time- step. On 
the other hand, injecting the energy before the hydro calcu- 
lation would make the current time-step inconsistent with 
the new hydrodynamical state of the particle. If the cur- 
rent time- step is overestimated, the integration could lead 
to unphysically too strong feedback effects. We also want 
to warn the reader against any feedback implementation 
that would inject the energy as a rate to be applied over 
a given length of time. Modifying either the rate of entropy 
change (thermal fedback) or the hydrodynamical accelera- 
tion (kinetic fedback) would not conserve the input energy 
if the kick operator (from the leap-frog integration scheme) 
is applied on two successive steps of different length. 
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Name 


Scheme 




/step 




^[h,k] 


Ener^ 
Thermal 


(%) 

Kinetic 


Time (: 
Thermal 


min) 
Kinetic 


SF,Dnv-n-D 

ij ejljkj V —yj(—LJ 


vjr lU L J cl 1 - 


\J \Jicx\jx5 






32 


0.93 


1.35 


139 


101 


SEDOV-L2-NU 


Limiter - 


No Update 


2 


0.0025 


32 


4.0 X lO'^ 


43.46 


119 


35.4 


SEDOV-L2-U 
SEDOV-L4-U 

SEDOV-L8-U 


Limiter - 
Limiter 

Limiter - 


Update 
- Update 

Update 


2 
4 

8 


0.0025 
0.0025 

0.0025 


32 
32 

32 


2.37 
3.03 

4.85 


2.19 
2.86 

5.28 


24.6 
21.8 

21.0 


20.8 
18.2 

17.7 


SEDOV-L4-U-r; 
SEDOV-L4-U-N4 


Limiter - 
Limiter - 


Update 
Update 


4 
4 


0.025 
0.0025 


32 

4 


3.90 
4.12 


4.28 
4.33 


21.0 
26.4 


16.8 
25.7 



Table 2. Energy conservation estimate and computational time for Sedov's blast wave tests with a limiter time-step scheme are compared 
to a reference simulation using global time-steps, when the injection of energy is delayed (D). The reference simulations, including our 
recommended set of parameters when using the limiter, are highlighted in bold characters. Names are given to each simulations according 
to the integration scheme used and when numerical parameters differ from the run of reference. Energy conservation estimates (in percent) 
are given at the end of the runs, at time t = 0.04 after energy injection. The wall-clock time (in minutes) is for runs on 8 processors. 



Considering the problems mentioned above, we recom- 
mend injecting instantaneously the feedback energy just be- 
fore the computation of the next time- step, whilst, taking 
into account the following two actions to preserve the accu- 
racy of the time integration. Firstly, to ensure the fast infor- 
mation transfer of the change of energy in the medium, we 
use the time- step limiter proposed by SM09. In Appendix C 
we describe an implementation of this limiter in GADGET-2 
which preserve the time- step synchronisation. Secondly, to 
ensure that the computation of the time-step following the 
explosion takes properly into account the local hydrodynam- 
ical change, we impose that the particles which receive the 
energy update their signal velocity and hydrodynamical ac- 
celeration. Doing so, we impose for these particles an update 
of the computation of their next time-step that will lead 
to an update of their hydrodynamical properties right af- 
ter the explosion time (see Appendix D for the detail about 
this time-step update). We will now show that these actions 
are essential to be considered altogether in order to preserve 
the concordance of feedback methods when using individual 
time-steps. 

3.2 Sedov's blast wave test 

In this section we describe the Sedov's blast wave tests per- 
formed with an individual time- step scheme. Before showing 
the results of our simulations, we mention the differences 
with the setup used in Sec. 2.1. 

3.2.1 Simulation setup 

Starting with the same glass-like uniform conditions pre- 
sented in Sec. 2.1, all particles are evolved over global, back- 
ground time-steps that are of the order of Atback ~ lO""^. 
Since all particles are synchronised from the beginning of 
the simulation, the issues about an explosion occurring in 
the middle of active steps or about neighbouring particles 
being initially on different time-bin levels, will not be ad- 
dressed in this section. We refer to Sec. 3.3 for an analysis 
of these effects on the long-term evolution of the medium. 

To avoid a pre-defined population of the time-bin lev- 
els in these Sedov's test simulations, the injection of a total 



energy of = 1 (either in thermal or kinetic form), is de- 
layed by a few of the background steps. In order to focus on 
feedback processes similar to SN explosions or BH activity, 
we have set the total initial internal energy of the system to 
^0 = 1- Given the resolution N — 128^ of the simulations, 
the initial energy contrast between the heated/kicked parti- 
cles and the background ones is of the order of ~ 10^. All 
tests presented in this section have been run until t — 0.04 
after the explosion time. 

It is interesting to remind here that the amount of in- 
jected energy constrains the time-steps that follow the ex- 
plosion. Since the signal velocity (see Monaghan 1997) at 
the explosion location is related to the energy injection as 
Vsig oc ^Ju^^ we can already anticipate through the Courant 
criterion that, for the energy contrast considered here, the 
step of the heated/kicked particles will be ~ 10"^ times 
smaller than the background step. This will correspond to an 
abrupt drop of about 10 levels in the hierarchy of time-bins. 

To analyse the behaviour of our time- step scheme, we 
use a set of simulations (see Table 2) where combinations 
of the integration techniques and numerical parameters are 
investigated. Here we use again a simulation with a global 
time-stepping scheme as reference. Then, we test the limiter 
technique without the time-step update. Finally, we enforce 
the time-step update, as describe in Appendix D. With the 
latter setup, we estimate the impact of the limiter parame- 
ter /step , as well as the time integration efficiency parameter 
Tf (from Eq. B3), for both the thermal and kinetic feedback 
methods. With the aim of studying two different initial en- 
ergy distributions, we also compare simulations where the 
energy is injected over a different number of particles. 

In order to quantify the accuracy of the different meth- 
ods, we estimate the input energy conservation error at 
t = 0.04 for each simulation using Eq. 3. To estimate the 
performances of the different time- step schemes, we also list 
the total running time of each test in Table 2. 

3.2.2 Simulation results 

We present in Fig. 5 the early state of the medium at 
t = 0.004, for both thermal (left column) and kinetic (right 
column) energy injection, when using the time-step limiter 
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Figure 5. Sedov's explosion test for thermal (left column) and ki- 
netic (right column) energy injection, when the limiter is applied 
without the time-step update. Results are given soon after the 
explosion at time t = 0.004. We see that for the thermal imple- 
mentation of feedback a stable shell has developed far ahead from 
the similarity solution because the large violation of energy con- 
servation largely increases the system total energy (see Table 2). 
In the kinetic feedback case inter-particle crossing is evident. No 
density contrast has yet developed given that kicked particles have 
traveled far from the explosion site. These results illustrate how 
the lack of a prompt response of the medium leads to two distinct 
effects from the two feedback schemes. Time animation available 
as supporting online material. 



without any update at the explosion time. This corresponds 
to the time integration that has been illustrated in the ap- 
pendices in Fig. 11-b). For those runs, we have chosen the 
conservative value of /step = 2 in order to ensure the fastest 
information propagation between neighbouring particles. It 
is striking to see how the two feedback methods give ex- 
tremely deviant but very different results. In the thermal 
case, heated particles are found close to the analytical blast 
radius but since the first step after the explosion is much 
too long, the acceleration given to the surrounding medium 
produces an important excess of kinetic energy. Eventually, 
the propagation of this excess of energy will lead to a viola- 
tion of energy conservation to a level of ^ 40, 000 % at the 
end of the run. 

For the kinetic injection of energy, we see that the 
isotropic nature of the blast is destroyed. Indeed, the lack 
of an immediate response of the medium makes kicked- 
particles cross their neighbours before they partially ther- 
malise their momentum. Later on, they will eventually start 
to expand individual bubbles around them. Even with the 
use of the limiter, some particles are able to travel up to the 
border of the volume before interacting with the medium. 
Here, the energy creation reaches a level of ^ 43%. This two 
pictures clearly show the importance of properly computing 
the time- step following the energy injection event. 

In Fig. 6, we show the results for the simulations in 



Figure 6. Sedov's explosion test for thermal (left column) and 
kinetic (right column) energy injection, when both the limiter and 
the time- step update are applied. Compared with the previous 
figure, results are shown at later time t = 0.02. We show here that, 
for our recommended set of parameters, the two feedback schemes 
are able to reproduce precisely the similarity solution and hence 
provide again concordant results as in Fig. 2. This demonstrate 
that both a fast information transport and a prompt response to 
the explosion are needed in order to resolve strong feedback events 
with individual time-step integration. Time animation available 
as supporting online material. 



which we additionally enforce the time- step update. Here we 
use the fiducial value of /step = 4. The shell position and the 
radial density profile are in extremely good agreement with 
the analytic solution. Sedov's solution is remarkably well 
represented for both feedback methods, reproducing again 
the concordance with a conservation of input energy close 
to 3 % at the end of the runs. 

Regarding the parameter study shown in Table 2, we 
can see that increasing the value of /step slowly enhances 
the energy violation, while keeping the computational time 
nearly constant. Using a less conservative accuracy param- 
eter r] = 0.025 increases also the energy creation but both 
thermal and kinetic simulations are completed faster. We 
note that even distributing the energy over a smaller num- 
ber of neighbours, which provides a larger energy contrast, 
gives an acceptable energy conservation. 

In Fig. 7, the time evolution of the energy conservation 
for the runs using the fiducial value of /step = 4 is compared 
with the reference simulations using the global time-step 
scheme. We notice that the ratio of energy increase between 
matching models is fairly constant. This behaviour tells us 
that the convergence between the two feedback methods is 
achieved since the very beginning of the expansion of the 
blast. Moreover, we see in this figure that all simulations 
experience first a sharp increase of the energy creation be- 
fore reaching a stable level of energy conservation. Indeed, 
the rate of energy violation is high right after the explosion 
(when the energy needs to be absorbed by the surround- 
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Name 
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Energy (%) 
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4 

8 
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0.0025 
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32 
32 

32 


2.01 
2.16 

2.27 


1.27 
2.60 

3.26 


3.63 
3.41 

3.35 


3.60 
3.39 

3.34 


HALO-L4-U-r; 
HALO-L4-U-N1 


Limiter - Update 
Limiter - Update 


4 
4 


0.025 
0.0025 


32 
1 


4.05 
2.60 


5.51 
0.95 


2.40 
3.44 


2.38 
3.47 



Table 3. Energy conservation estimate and computational time for the off-centre explosion tests, where the energy injection is delayed 
(D) after the start of the simulation. The reference simulations, including our recommended set of parameters when using the limiter, 
are highlighted in bold characters. Names are given to each simulations according to the integration scheme used and when numerical 
parameters differ from the run of reference. Energy conservation estimates (in percent) are given at the end of the runs, at time t = 0.04 
after energy injection. The wall-clock time (in hours) is for runs on 8 processors. 
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Figure 7. Time evolution of the energy conservation for the most 
relevant Sedov's blast wave tests. Solid lines show the results for 
the thermal injection of energy. Dashed lines represent the evolu- 
tion of the kinetic feedback runs. 



ing of the injection region) but as soon as the blast be- 
comes stable and starts to expand, the propagation of the 
error becomes very small. This happens approximately for 
t > 0.004. We also confirm from this analysis that the time 
accuracy parameter 77 needs to be chosen small enough in 
order to get the best convergence with the feedback meth- 
ods. Finally, regarding the spatial energy distribution, it is 
clear that concentrating the injection over a smaller num- 
ber of particles gives a higher energy jump (meaning a larger 
time- level drop), and hence produces a slightly worse energy 
conservation. 

We conclude that, when applying the proposed inte- 
gration scheme, there are no dramatic differences either in 
concordance, or in energy conservation accuracy, or in per- 
formances for the set of parameters we explored. 



3.3 A more realistic test 

In this section we simulate an explosion event in a self- 
gravitating gas halo. The initial conditions are similar to 
the ones in SM09. However, the injected energy is placed off- 
centre to follow the blast evolution in the presence of pres- 
sure and density gradients. In cosmological simulations of 
structure formation, the sources of mechanical and thermal 
feedback (SNe and BHs) are generally situated in a similar 
environment. Though idealised, the setup gives a qualitative 
and quantitative view of the above scenario. 



3.3.1 Simulation setup 

We create a spherical particle distribution of density profile 
p (X by spatially remapping a uniform, spherical distri- 
bution of particles through the radial transformation 

ro\d,i rnew,2 = (~^^) ^' (^^^ 

where R = 1 in the same system of units used for Sedov's 
blast wave problem. The sphere is cut out from a uniform 
distribution of A/" = 128"^ particles in a cubic volume of side 
L = 2. The total mass of the gas sphere is unity. We assign an 
internal energy of 0.05. The system is evolved including self- 
gravity up to time t — 3 when relaxation has already taken 
place (Evrard 1988). We chose the gravitational softening 
e — 3.125 X 10""^ and the artificial viscosity parameter a — 1. 

We start the simulation tests using the relaxed gas halo 
as initial conditions, and inject either thermal or kinetic en- 
ergy in a region centred on [0.05, 0, 0]. The off-centre explo- 
sion is generated by enhancing the energy of particles by 
U[h,k] as described in Sec. 2.1. We avoid the first simulation 
step in which all particles are active, and delay the energy 
injection to the time t^ — 10~^. For the sake of simplicity, 
in the whole analysis, we refer to the explosion time as the 
initial time to = t^. Any other time is relative to to- To re- 
move another bias, we also ensure that at injection time all 
particles receiving energy are inactive. Doing so, we can ap- 
ply our time-stepping scheme to the most generic situation, 
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Figure 8. The off-centre explosion in a self- gravitating gas sphere 
for thermal (left column) and kinetic (right column) energy in- 
jection. In all pictures we show the gas density at time t = 0.02 
in a slice of thickness Az = 0.1 centred on the box origin. The 
white dots mark the position (at the same time) of all particles 
that received thermal or kinetic energy. Top row: only the time- 
step limiter is applied using the conservative value of /step = 2. 
The energy violation is severe in the thermal case and the bubble 
have blown away a large fraction of the gas halo. Bottom row: 
only the time-step update is enforced. The behaviour is similarly 
wrong in both cases. The halo atmosphere is disrupted and no 
expanding bubble forms. Time animation available as supporting 
online material. 

but more importantly, put it also under the least favourable 
conditions. 

The mass-weighted background internal energy at the 
position of the explosion is u = 1.75, whereas the mass- 
weighted velocity is = 4.43 X IQ-^ The average internal 
energy and velocity are measured in a spherical shell of thick- 
ness 0.004 and average radius 0.05. In the case of a single 
heated particle, the jump in energy would be :^ 5 x 10^ times 
its initial value. If only one particle is kicked, its momentum 
change would be ^ 3 x 10^ times the background initial 
value. 

We compute the energy conservation as the relative 
variation of the total injected energy as in Sec. 2, and includ- 
ing the potential energy. The simulation parameters together 
with the energy conservation estimate and the running time 
are listed in Table 3 where the reference simulations are 
highlighted in bold. 

3.3.2 Simulation results 

We show in Fig. 8 the off-centre explosion in a self- 
gravitating gas sphere for thermal (left column) and kinetic 
(right column) energy injection. In all panels we show the 
projected gas density at time t = 0.02 in a slice of thickness 
Az = 0.1 centred on the box centre. The white dots mark 
the position (at the same time) of all particles that initially 
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Figure 9. The off-centre explosion in a self- gravitating gas sphere 
for thermal and kinetic energy injection when time-step limiter 
and update are applied. The plots are as in Fig. 8. Top row: 
the energy is injected on 32 particles in thermal (left panel) and 
kinetic (right panel) form. Impacted particles are made active and 
the time-step of surrounding particles has been corrected. The 
results are qualitatively identical, showing concordance of the two 
feedback methods. Bottom row: all available energy is injected 
into one particle. The results excellently agree with each other and 
with the cases in which the energy is injected into 32 particles. 
Time animation available as supporting online material. 

received the energy. We first note the effect of the pressure 
and density gradients. In all pictures, the heated/kicked par- 
ticles have moved from the location of energy injection along 
the direction of decreasing gradients (x-axis). 

In the upper panels, we see the result of correcting the 
feedback scheme by making impacted particles active. Al- 
though heated/kicked particles are integrated accurately, 
their neighbours do not react soon enough. Indeed, be- 
fore neighbours are active again, heated particles accelerate 
along the decreasing pressure gradients, and kicked parti- 
cles build thermal energy through shocks. Once the neigh- 
bours feel the energy perturbation, they receive large accel- 
erations. Integrating over a too long time- step gives them 
unrealistic momentum values, and in a few steps they are 
ejected from the vicinity of the explosion. Without the use 
of the limiter, the propagation of the effect from neighbour 
to neighbour builds up energy to such levels that its con- 
servation is violated by - 57, 000 % and - 170, 000 % for 
thermal and kinetic energy injection, respectively. It is in- 
teresting to see that both feedback implementations lead to 
the complete disruption of the halo atmosphere, which lost 
its initial spherical shape. This confirms the results of SM09 
about the problems encounter with the standard individual 
time-stepping scheme and shows the importance of ensur- 
ing a smooth transport of the information from particle to 
particle. 

In the lower panels of Fig. 8, we show the tests using the 
time-step limiter but without the time-step update. Though 
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the limiter factor is set to the conservative value of two, the 
energy conservation is still largely violated. 

With thermal feedback (bottom- left panel), we recover 
the behaviour discussed in Sec. 3.2.1. The blast wave ex- 
pands to large radii because thermal energy is overproduced 
within the hot bubble. Indeed, energy conservation is vio- 
lated by ~ 3, 000%, and the total thermal energy of the sys- 
tem is boosted to a factor of 20 or more the input one. At 
the final time, half of the thermal energy has been converted 
into momentum, and the total kinetic energy is two thirds of 
the total. Simulations of strong BH and SN feedback should 
be carefully checked against this behaviour. 

If kinetic energy is injected (bottom-right panel), the 
picture is similar to that discussed in Sec. 3.2.1. Inter- 
particle crossing is clearly recognisable in the picture. The 
particles with highest velocities (the closest to the injection 
position) have crossed several smoothing lengths and shock- 
heated the gas farther away. We note that within the blast 
wave several small bubbles are created, and dense gas is 
entrained, possibly decreasing the shell mass. Though less 
severe than the thermal case, energy conservation is vio- 
lated by ^ 20%, which, in cosmological simulations may 
lead to overestimating the effect of SN winds. In summary, 
we showed that neither the update nor the limiter of im- 
pacted particles time-step are able by themselves to conserve 
energy within a reasonable accuracy level. 

We show in Fig. 9 the off-centre explosion test when 
both the time-step limiter and update are applied (all plots 
are as in Fig. 8). We show in the upper panels the tests with 
^[h,k] — 32, and in the lower ones the test with n[h,k] = 1- 
The agreement between the same model with different n[h,k] , 
and between different models, is striking. The evolution 
of the hot bubble is very similar in all cases, expanding 
and buoyantly moving from the place of energy injection 
to larger radii. The reference simulations with the global 
time- stepping integration show a conservation error < 2%, 
whereas the reference simulations with individual time-steps 
(/step = 4, n[h,k] = 32) have conservation error < 3%. We 
warn the reader that the gravitational force calculation is 
accurate at ^ 1% level and could contribute to the errors 
listed in Table 3. 

Nonetheless, the lower panels of Fig. 9 show that con- 
cordance of thermal and kinetic feedback is achieved even 
in the case of heating/kicking just one particle with all the 
available energy. 

The time evolution of the energy conservation relative 
error is plotted in Fig. 10. We plot some of the most rel- 
evant thermal (solid lines) and kinetic (dashed lines) feed- 
back tests. Each line colour refers to simulations that were 
run with the same numerical parameters. The agreement be- 
tween kinetic and thermal feedback is excellent for the global 
time-stepping scheme (black lines), proving that a prompt 
response of the system to the energy perturbation is neces- 
sary to achieve concordance of the two methods. Applying 
the proposed time-step scheme (red lines, reference model), 
the integration slightly looses accuracy, but simulations get 
a considerable speed up by a factor of ^ 8 (see Table 3). 
The agreement between the two feedback schemes is again 
excellent. 

We show in blue the tests with 77 = 0.025. Given the 
lower time integration accuracy, energy violation naturally 
reaches higher values, but the ratio of relative errors is simi- 
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Figure 10. Time evolution of the energy conservation for the 
most relevant tests of the explosion in a gas halo. Solid lines 
show the results for the thermal injection of energy. Dashed lines 
represent the evolution of the kinetic feedback runs. 

lar to that of the reference tests. We see that the dashed blue 
curve goes to higher values right after the explosion. That 
shows that it is crucial to properly capture the conversion of 
momentum through viscous forces at the very earliest time. 
At the final time, the relative error is 4 and 5.5% for thermal 
and kinetic schemes, respectively. 

Finally, simulations where only one particle is 
heated/kicked (green lines) also show a good concordance 
of feedback methods. However, injecting kinetic energy leads 
this time to a loss of total energy. This can be explained be- 
cause all particles close to the explosion are inactive when 
energy is injected. Since all the energy is concentrated, the 
kicked particle has the time to move for a few steps while 
decelerated by viscous forces and before the surrounding 
medium reacts. During this time, there is no energy transfer 
to neighbours which lead to the loss of energy. This problem 
could be solved by making all neighbours active at the same 
time of the impacted particles. We have tested this hypothe- 
sis, and obtained again the concordance. In any case, energy 
is still conserved with good accuracy. 



4 CONCLUSIONS 

In this work we investigated the concordance of thermal and 
kinetic feedback methods in the absence of cooling. We nu- 
merically demonstrated that the two methods are equiva- 
lent using simulations of Sedov's explosion with a global 
time- step integration scheme. We argued that this result is 
expected from the hypothesis of Sedov's problem: a large 
amount of energy is instantaneously injected in a small vol- 
ume. The blast physical properties are given by the integra- 
tion of hydrodynamics conservation laws which also give the 
fraction of energy in thermal (72%) and kinetic (28%) form. 
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We focused on the conservation of input energy which 
we achieved within 1% in the reference simulation. In test 
simulations with varying input energy, the energy conserva- 
tion relative error flattens with time to roughly the same 
value, and is independent of the energy for a constant arti- 
ficial viscosity parameter. Decreasing the artificial viscosity 
parameter results in a larger conservation error. We also in- 
vestigated the evolution of the thermal and kinetic fractions 
of the input energy. We found excellent agreement with the 
expected values derived from Sedov's analytic solution. 

Using the global time-step scheme in cosmological simu- 
lations is computationally expensive. To reduce the compu- 
tational costs most of the numerical codes use an adaptive, 
individual time-step scheme. SM09 showed that to achieve 
good energy conservation of feedback energy a time-step lim- 
iter is necessary. We implemented a similar limiter in gad- 
get- 2 taking into account time- step synchronisation and 
hierarchy, and performed simulations of the Sedov's blast 
wave. We showed that the limiter does not give accurate 
energy conservation if the explosion energy is injected at 
random time and a maximum value of the time-step for am- 
bient gas particles is not enforced. 

We proposed a solution to the problem: the system must 
not only propagate the information rapidly but should also 
promptly react to the energy input. That is done by updat- 
ing the time-step of particles receiving energy accordingly 
to their new hydrodynamical state. 

We tested the modified time integration scheme on the 
Sedov's blast wave, and obtained good energy conservation 
accuracy and concordance of methods. We applied both the 
limiter and the time-step update to a more realistic case of 
an explosion in a self- gravitating gas halo, and showed that 
concordance and accuracy are achieved also in the extreme 
case where all the available energy injected into one particle. 

In the test simulations presented here, we did not con- 
sider radiative cooling, but we are aware of the potential 
breaking of concordance when such processes become im- 
portant. Cooling being a non- linear physical process, the 
amount of radiated energy before reaching energy partition 
will depend on the specific feedback scheme and the conse- 
quent gas temperature evolution. In the thermal feedback 
case, the bubble temperature decreases by adiabatic expan- 
sion and conversion of thermal energy into momentum. In 
the kinetic case, the temperature increases through shocks. 
It is thus important to reach the same energy configura- 
tion before any significant cooling loss. As mentioned by 
Dalla Vecchia & Schaye (2008), if the gas is heated to a tem- 
perature where the cooling time becomes longer than the lo- 
cal dynamical time, radiative over-cooling may be prevented. 
These considerations will be addressed in more detail in an 
upcoming work. 

We summarize this work as follows: 

• concordance of thermal and feedback methods arises by 
accurate time integration and is expected from theoretical 
arguments; 

• keeping all other numerical parameters fixed, the max- 
imum energy conservation relative error is roughly constant 
when varying the input energy, and at the same level with 
both feedback schemes; 



• high artificial viscosity coefficient enables a fast conver- 
sion of thermal energy to momentum, and hence permit to 
numerically converge to a stable solution in a shorter time; 

• concordance and energy conservation can be achieved 
with a time-step limiter, provided that the system promptly 
reacts to the input of energy; 

• in cosmological simulations with strong kinetic and/or 
thermal feedback from SNe and BHs, one should take into 
account accurate energy conservation before inferring any 
results dependent on feedback processes. 
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In the following appendices we describe all the relevent 
part of the time integration implementation of GADGET-2. 
The reader who is not familiar with the leapfrog integra- 
tion may read Appendix A. There, we also remind about 
the individual time-step hierarchy on which the leapfrog is 
applied. In Appendix B we review the two criteria used to 
define the individual time-steps and explain the modification 
we made to match our need of accuracy. In Appendix C, we 
proceed to the description of the implementation of the time- 
step limiter in GADGET-2. Finally, we discuss in Appendix D 
our time-step update which enables to achieve a prompt re- 
sponse of the system to any type of energy perturbation. 



that characterise any particle as follows: 
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In the above equations, Xn, Vn, an, Pn, Pn, Sn and Sn are 
respectively, position, velocity, acceleration, pressure, den- 
sity, entropic function and its rate of variation at time tn- 
It is interesting to note here that the acceleration and the 
rate of entropy change need to be computed before the K 
operator is applied (which implies an update of the system), 
while the drifting is a free-motion operation that uses the 
velocity and entropy values from the previous kick. 

In cosmological simulations of galaxy formation, where 
a large range of time scales is required, the use of an in- 
dividual time-step scheme enable to reduce the number of 
force computations per step to only those needed by a frac- 
tion of the system. The KDK integrator is then used over 
small time-steps, 6t, for those particles, while the other par- 
ticles experience fewer kicks over longer steps. However, 
given the pair- wise nature of the forces, the potential is re- 
quired to evolve on the shorter time scale, in order to pre- 
serve the accuracy of the force computation. As long as the 
larger time- steps can be subdivided into several shorter ones. 
At = J2k ^^k, the splitting of the drift operator for the long 
time-step particles ensure that the positions of all particles 
are updated for any force calculation: 



(A5) 



APPENDIX A: LEAPFROG INTEGRATOR 

For a given time interval At = tn+i — tn, the KDK leapfrog 
integrator used in GADGET-2 (first presented in this form by 
Quinn et al. 1997) is defined by the following time evolution 
operator: 



U{At)=K'(^)D{At) K"^^* 



where the kick (K) and drift (D) operators are used to 
evolve both the dynamical and hydrodynamical quantities 



If the individual time-step scheme allows to integrate 
each particle according to its local dynamical state, one side 
effect is that the positions of all the particles of the system 
need to be drifted in order to compute the forces at a given 
time. To optimise the efficiency of the leapfrog integrator, 
one has to impose that several particles share exactly the 
same simulation time. This can be achieved by discretising 
the time-steps, as proposed by Hernquist & Katz (1989) and 
Makino (1991), in a hierarchy of power of two subdivisions 
of the global simulation time, tgiob: 

tglob 

where /c > is an integer that define the level of the corre- 
sponding time-bin. Particles are allowed to move to a smaller 
time- step, but can only jump up to the next larger time- step 
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a) Hierarchical Individual Time-steps 
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b ) Time-step Limiter 
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c ) Time-Step Limiter + Time-step Update 
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Figure 11. Timelines for a synchronous time integration scheme. Energy is given to particle i, and particle j is one of its neighbours. 
Solid lines represent time intervals where the state of particles is re-computed, while long-dashed lines correspond to "fictitious" steps. 
Circles represents the "active" times of particles. Labels given on the right side of each sketch show the levels of the time-bin hierarchy, 
a) Individual time-steps: particles are only able to adapt their time-step when the simulation time equals their end-step time (when they 
become active). Therefore, heated/kicked particles may complete a significant number of steps before their neighbours become active as 
well, b) Time-step limiter: neighbouring particles communicate to each other the length of their time-steps, and keep the ratio of long to 
short time-steps not larger than a factor of four (/step = 4). However, particles still need to be active before their neighbours adjust their 
time-step accordingly to the limiter criterion, c) Time-step limiter and time-step update: in the scheme introduced in this work, we make 
sure that heated/kicked particles become active at the time of energy injection, and adjust their time-step accordingly to the amount of 
energy they receive. When applied in combination with time-step limiter, impacted particles and their neighbours can promptly react to 
the change of energy in the medium. 



every other step (see panel a of Fig. 11 for an illustration). 
An additional advantage of the hierarchical scheme is that 
the previous constraint leads to the synchronisation between 
successive time levels which helps to balance the distribution 
of particles into time-bins. 

In GADGET-2, the length of the particle time-step is 
given by two variables storing the beginning, tbeg, and end- 
ing, tend, time of the current particle step. All particles with 
tend equals the current simulation time are considered ac- 
tive. This means that their dynamical and hydrodynamical 
states are updated by the following operator: 

f , grav I hydro 

^(tend) : 1 1"^^ 5^"' " ^"'^ + ^""^ (A7) 

I ^^beg ^ J^end 

where the acceleration, a, describes the motion of the 
fluid element, and takes also into account the artificial vis- 
cous forces generated by the shocks that develops into the 
medium. As a consequence of the local viscosity, the mo- 
mentum of the gas is converted into entropy at a rate given 
by S. However, before applying the operator A, some other 
quantities need to be updated, such as the smoothing length, 
and the local density and pressure. The reader should refer 
to the original GADGET-2 paper for more details about the 
specific expression of these quantities. 

To conclude this reminder of the time integration 
scheme used in GADGET-2, we want to remark that the two 
kick operators and K'' given in expression A2 and A4, 
will use the acceleration and the entropy rate computed at 
the same time, tend, when applied in two successive steps. 
The implementation of the leapfrog method in GADGET-2 
makes use of this property to save some computational op- 
erations. At the end of each active step, the length of the 
next time-step, Atnext, is computed before applying the kick 
operator from the current half- step to the next half- step. 



Therefore, the time evolution operator can be rewritten as: 

[/ = D(tbeg ^ tend) (A8) 
^(tend) K ^tbeg H > tend H ' 

where D is applied at any simulation step that sub-divide 
Atnow, while both A and K are only evaluated at the end 
of the individual steps for active particles. 



APPENDIX B: TIME-STEP CRITERIA 

Time-steps are chosen as the smallest given by the Courant 
criterion and the acceleration criterion. The Courant time- 
step criterion is defined by the signal velocity method, t'sig.i, 
(Monaghan 1997) which estimates the speed at which the 
information propagates at the spatial resolution scale, 

Atc,i = C^. (Bl) 
where C is the Courant factor and hi the smoothing length.^ 

The acceleration time-step criterion is defined in GAD- 
GET-2 as 

Ataccz = aA — , (B2) 

V Cii 

where 77 is the accuracy parameter, e is the gravitational 
softening, and ai is the particle acceleration.^ Though this 

^ Note that in GADGET-2 the smoothing length h is the size of the 
SPH kernel, whereas other authors define the kernel of size 2h. 
Hence, the fiducial value of the Courant factor used in GADGET is 
half of the one found in other SPH codes. 

^ In pure hydrodynamic simulations, the acceleration is given 
only by hydrodynamical forces while, in simulations that also 
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criterion gives accurate integration in most applications, it 
is crucial that the initial, large hydrodynamical accelera- 
tions, encountered under strong energy perturbations, are 
properly captured. However, since we allow the individual 
smoothing length hi to be smaller than the gravitational 
softening, the above criterion could overestimate the length 
of the time-step. Following Monaghan (1997), we thus mod- 
ified the acceleration criterion for SPH particles as follows: 



2 mm{hi, e) 



(B3) 



In addition, as reported by Wuyts et al. (2010), the 
value of the accuracy coefficient, 77, may play a role in cap- 
turing strong shocks, and hence in the convergence of the re- 
sults. This motivated the accuracy tests presented in Sec. 3.2 
and 3.3 which made us choose 77 = 0.0025 as fiducial param- 
eter for this work. 



APPENDIX C: TIME-STEP LIMITER 

We generalise the implementation of the time-step limiter 
in order to deal with energy injection happening at random 
time and affecting active and/or inactive particles. We de- 
scribe here the implementation of the time-step limiter into 
GADGET-2. The main idea behind the limiter is to ensure 
that every neighbouring particle j is integrated over a time- 
step that is shorter than or equal to a multiple of particle i 
time- step: 



Atj < /stepAti 



(Cl) 



In the context of the kick-drift-kick (hereafter, KDK) 
leapfrog integration presented in Appendix A, particle i, 
which experiences a change of its energy content in the 
middle of its current time- step, will only be able to re- 
act to the explosion when it becomes active. This is il- 
lustrated in Fig. 11 a). At this time, and assuming that 
the energy change is carefully taken into account (see dis- 
cussion in the next Appendix), the particle next time-step, 
Ati, will be brought to a finer level allowing a more accu- 
rate computation of the hydrodynamics. Later on, when the 
medium permits it, the particle time-step will eventually in- 
crease according to the synchronisation scheme. As shown 
in Fig. 11 a), the neighbouring particle j will have to wait 
a long time before becoming active and acknowledging the 
energy change. Therefore, to make sure that all the neigh- 
bouring particles will react to the explosion, particle i has 
to communicate them its next time- step. 

The implementation of the limiter is illustrated in 
Fig. 11 b) and proceeds as follows. We consider particle i 
which receives energy at time , and particle j which is one 
of its neighbours. When active, particle i informs neighbour 
j of its next time- step. If the neighbouring particle j is active 
and the inequality in Eq. Cl is not satisfied, its time-step 
will be shortened to the value Atj — /stepAt^. If particle j 



is inactive, the time-step limiter criterion will ensure that it 
will be recomputed no later than: 



+ /stepAti , 



(C2) 



where tend,i is the end-time of particle i and Ati is its next 
time-step.^ If particle j time-step ends later than the limiter 
time, tend,j > tiim, we change it to the new value: 

tend,j — ^beg,j + kf step Ati , (^3) 

where k is the largest integer that verifies the condition 



theg,j ~\~ k f step Ati ^ ^lin 



(C4) 



This ensures that particle j is computed within tum- A new 
time- step, Atj = f step Ati, is thus assigned to particle j and 



the beginning of the step becomes t^g 



end,_7 



At',. 



This is done to preserve the synchronisation in the time-step 
hierarchy. Finally, because both the beginning and end-step 
of particle j have been modified, velocity and entropy must 
be made consistent with the new time- step, by interpolating 
their values from the middle of the original time-step to the 
middle of the new one. This can be done by applying the 
following kick operator (see Appendix A) 



K tbeg,.7 + 



At,- 



+ 



At', 



(C5) 



Since the particle is inactive, the acceleration and the rate 
of entropy change used in the interpolation are the ones 
computed at the beginning of the previous active step. 

Regardless of the implementation of the limiter, we can 
clearly see in Fig. lib) that, if the explosion occurs in the 
middle of the particle i time- step, a considerable amount 
of time passes before both particles are recomputed accord- 
ingly to this event. In order to prevent this undesirable be- 
haviour, we present in the next section our suggestion to 
correct the timeline of impacted particles. 



APPENDIX D: TIME-STEP UPDATE 

The basic idea behind our time-step update can be divided 
into two specific actions: i) to make sure that the computa- 
tion of the step following the explosion takes into account 
the local change of the energetics; ii) to ensure that impacted 
particles react as soon as the energy is injected. 

In the following, we consider an instantaneous injection 
of energy, where the first time-step following the explosion 
event defines how far the information can travel before the 
medium starts to absorb the energy. It is thus crucial that 
the time-step of impacted particles does not violate either 
the Courant or the acceleration criterion. In our implemen- 
tation, we therefore update the maximum signal velocity and 
compute the acceleration of impacted particles just after the 
call of the feedback routine. The acceleration will be used 
only in the calculation of the time-step, and does not substi- 
tute the one computed at the beginning of the step. Given 
the pair-wise nature of the hydrodynamical forces, the time- 
steps of active particles in the neighbourhood of the energy 
injection area are conveniently updated at the same time. 



include self-gravity, the particle acceleration is the sum of the 
contributions from gravitational and hydrodynamical forces. 



^ The time tend,i is actually the current simulation time, since 
particle i is active. 
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The computational cost of an additional neighbour search 
is only a negligible fraction of that required for updating 
the densities and the hydrodynamical forces. Finally, this 
method has the advantage of naturally taking into account 
the effect of multiple sources of feedback. 

The second action that is needed for a prompt response 
of the medium, is to make active all the particles inside the 
explosion area. As shown in Fig. 11 c), this rather trivial op- 
eration allows particle i to shrink its time-step straight after 
the explosion time, t^. We can see from the figure that, when 
used jointly with the time-step limiter, both the impacted 
particles and their neighbours are now integrated on time 
scales that are able to correctly handle the upcoming evolu- 
tion of the blast wave. 

To briefly describe the implementation, we refer again 
to the KDK leapfrog scheme in Appendix A. At the time of 
explosion t^, the positions of all particles have already been 
drifted with the D operator. Since we want the particles to 
be recomputed just after the explosion, we have to apply 
again the kick operator to complete a "fictitious" step. The 
main difference with the kick correction used for the limiter 
being that the operator is here applied back in time: 

f theg,i + tend,i ^ ^beg,2 + ^end,2 

V 2 2 

where tend,i is the previous end-step time of particle i and 
^end,2 — is the current simulation time. 

Making the particle active at the explosion time allows 
particle i to recompute its next time- step Atnext,i according 
to the Courant criterion. Since this particle has not really 
been updated, the kick which closes this fictitious step, 
will use the acceleration and rate of entropy change com- 
puted at time theg,i- This is consistent with an instantaneous 
injection of energy, and prevents any loss of it. As illustrated 
in Fig. 11c), the heated/kicked particle is now able to handle 
its new energy content more accurately and can inform its 
neighbours immediately through the time-step limiter crite- 
rion. 




